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ABSTRACT 

To facilitate the study of black hole fueling, star formation, and feedback in galaxies, 
we outline a method for treating the radial forces on interstellar gas due to absorption 
of photons by dust grains. The method gives the correct behavior in all of the relevant 
limits (dominated by the central point source; dominated by the distributed isotropic 
source; optically thin; optically thick to UV/optical; optically thick to IR) and reason- 
ably interpolates between the limits when necessary. The method is explicitly energy 
conserving so that UV/optical photons that are absorbed are not lost, but are rather 
redistributed to the IR where they may scatter out of the galaxy. We implement the 
radiative transfer algorithm in a two-dimensional hydrodynamical code designed to 
study feedback processes in the context of early-type galaxies. We find that the dy- 
namics and final state of simulations are measurably but only moderately affected 
by radiative forces on dust, even when assumptions about the dust-to-gas ratio are 
varied from zero to a value appropriate for the Milky Way. In simulations with high 
gas densities designed to mimic ULIRGs with a star formation rate of several hundred 
solar masses per year, dust makes a more substantial contribution to the dynamics 
and outcome of the simulation. We find that, despite the large opacity of dust to UV 
radiation, the momentum input to the flow from radiation very rarely exceeds L/c 
due to two factors: the low opacity of dust to the re-radiated IR and the tendency 
for dust to be destroyed by sputtering in hot gas environments. We also develop a 
simplification of our radiative transfer algorithm that respects the essential physics 
but is much easier to implement and requires a fraction of the computational cost. 
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1 INTRODUCTION 

Supermassive black holes (SMBHs) are thought to inhabit 
the central regions of nearly all massive galaxies and have 
properties well-correlated with those of their host galaxies 



( Gebhardt et al.||2000[ [Ferra rese fc Merritt||2000[ [Tremaine 



et al. 2002 Novak et al. 2006; Giiltekin et al. 20091. The 



direction of the causal link between SMBH and host galaxy 
properties remains unclear, but the existence of such a cor- 
relation implies an intimate connection between the physics 
of SMBH growth and galaxy formation, spanning a factor 
of at least 10'* in length scale from ~ Ipc to ~ lOkpc. 

Quasars came to astronomers' attention due to their 
enormous electromagnetic output. Bright active galactic nu- 
clei (AGN) emit as much as 10''^ ergs"'^ over the electromag- 
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netic spectrum ( |Woo fc Uriy||2002| [Kollmeier et al.||2006[ ). 
A significant amount of careful work has been performed to 
study the physical effects of this extraordinary outpouring 
of electromagnetic energy and momentum on the material 
in the immediate vicinity of the AGN (e.g. Proga et al._2000[ 



2008[ [Kurosawa fc Proga|2009[ [Park fc Ricotti|2011[|20l"2 | 
As a result we are beginning to understand the radiation 
driven broad-line winds emerging from discs around the cen- 
tral SMBHs. However, relatively little work has been done 
regarding the radiative effects on the surrounding galaxy. 

There has been intense theoretical interest in the ability 
of AGN to affect the star formation histories and observed 



colors of galaxies (e.g. Croton et al.[[2006 Schawinski et al. 
20071. The relationships between star formation, and AGN 



activity, galaxy mergers, and secular processes have turned 
out to be complex and the parameter space is not yet fully 
mapped. What is clear is that gas added to the central re- 
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gions of galaxies will induce a central star burst and an AGN 
flare-up, but the triggering mechanism remains uncertain. 

Although the use of numerical simulations to study the 
interactions between black holes and galaxies has a long his- 
tory (e. g. [Binney fc Tabor|p95| [Ciotti fc Ostriker 1997|) 
work by pi Matteo et al.| (|2005[) and|Springel at al.('(|2005l 



set off an explosion of theoretical work. They studied the 
simplest reasonable fully specified model of black hole fuel- 
ing and AGN feedback in the context of a major mergers of 
gas-rich spiral galaxies and found that the AGN has a pro- 
found effect on the gas dynamics and star formation history 
in the merger remnant. 

During the combined AGN/starburst phase, analytical 



investigations by Thompson et al. ( 2005 1 and Murray et al 



1 2005 1 have shown the importance of radiative forces on 
dust. Subsequent numerical work ( [Debuhr et al.|2010)|2011[ 
2012 1 has verified that if the central ~ 100 pc of galaxies are 



assumed to have optical depths of 10 to 25 in the infrared, 
due to dust opacity, then radiative forces on dust have dra- 
matic dynamical consequences in the context of merging 
galaxies. 

As the investigations have become more detailed and 
accurate, it has become increasingly necessary to perform 
the radiative computations in a more rigorous fashion, care- 
fully determining both the energy and momentum input into 
the ambient gas. It is the purpose of this paper to outline 
some of the principles and techniques that may be used for 
this task. 

Our goal is to develop a numerical method to self- 
consistently solve the radiative transfer problem in galaxies 
in the presence of dust and radiation produced by stars and 
AGN. That is, we would like to solve for the specific inten- 
sity of the radiation field for a given distribution of radiation 
sources in the presence of scattering and absorbing media, 
and subsequently use the specific intensity to compute the 
forces on the scattering/ absorbing media. Studies of SMBH 
fueling and AGN feedback performed to date have either ig- 
nored radiative effects or else assumed a particular form of 
the photon field which is then used to compute forces on the 
gas. 

Although it would be ideal to obtain a perfect solu- 
tion for the radiation field in all cases, we will settle for a 
method that gives the correct asymptotic behavior in all rel- 
evant limits and interpolates reasonably between the limits. 
It should respect global conservation laws for energy and 
momentum. It must not be overwhelmingly demanding in 
terms of algorithmic complexity or computational resources. 

The paper is organized as follows: In Section[2]we sum- 
marize previous theoretical work on AGN feedback as well 
as some recent observational results. In Section [3] we dis- 
cuss the radiative transport equation in the context of the 
present problem. We describe the physics behind it, its most 
important properties, and the different approaches that can 
be used for the numerical solution. In Section |4] we describe 
the physics implemented in a two-dimensional hydrodynam- 
ical code designed to study black hole fueling and radiative 
feedback from AGN and star formation in the context of 
early-type galaxies. In Sectionlslwe present the main results 
of the simulations. Finally in Section [6] we summarize our 
conclusions. 



2 PREVIOUS WORK 

It has long been known that the bulk of the matter that 
enters black holes over the history of the universe does so via 
radiatively efficient accretion ( Soltan|[l982 Yu & Tremaine 
20021. The resulting photons are the clearest observational 



indication that SMBH accretion is taking place. 

There is a prodigious amount of energy and momentum 
available during Eddington-limited bursts consequent to ac- 
cretion, certainly enough to have a profound effect on the 
structure and dynamics of the interstellar medium. 



2.1 Mechanical feedback 

Many AGN are observed to have broad absorption lines in- 
dicating strong outfiows. Several physical processes oper- 
ating near the SMBH could plausibly launch these winds 
(e.g 



Proga et al. 20001 and they go under the heading of 



mechanical feedback regardless of the details of their origin 
( [Ostriker et al. 20101. They are thought to be somewhat 
collimated and are thus only observed in 1/4 to 1/2 of all 
bright AGN. Recently, quantitative estimates of the kinetic 



energy in these outflows have become available (Moe et al. 
Arav et al.||2012| and indicate that they contain 0.1- 



2009 



1% of the energy in the radiative component of the AGN 
output. 

This mechanical energy can couple to the ISM much 
more efficiently than the AGN photons, making this an at- 
tractive mechanism for AGN feedback. However, there are 
significant observational uncertainties in these kinetic wind 
energy estimates, the studies have only been carried out for 
a few objects, and there is no global constraint analogous to 



the Soltan ( 1982 \ argument that allows an estimate of the 



global mean efficiency of energy conversion to kinetic winds 



in AGN. If the Moe et al.|(2009J estimate of the kinetic lumi- 
nosity is correct and we assume all of the mechanical energy 
and momentum is effectively deposited into the ISM, then 
mechanical feedback will dominate over radiative feedback 
for column densities less than 10~^^ to 10~^^gcm~^. 



2.1.1 Simulations of mechanical feedback 

Many groups have carried out substantially similar studies 
of AGN feedback in the context of galaxy mergers where the 
black hole accretion rate is given by the smaller of the esti- 
mated Bondi rate and the Eddington rate, while the AGN 
feedback is assumed to be purely thermal energy injected in 
the vicinity of the black hole. This is the basic model devel- 
oped by pi Matteo et all (|2005l) andlSpringel et al.|(|2005l 



and applied extensively to compare with observations of the 
quasar luminosity function, galaxy colors, and more (e.g. 
Hopkins et"arf2 005 200 6l|2008[). Th e same basic model was 
also used by Johansson et al. (20091, and all authors found 



that AGN feedback has a dramatic effect on galaxies pro- 
vided that 5% of the bolometric AGN output is transformed 
into thermal energy within 100 pc of the black hole by me- 
chanical winds operating on smaller scales. 

However, observations of gas dynamics near AGN have 
indicated that the actual efficiencies are a factor of 5 to 50 
smaller than the value usually used in many of these simu- 
lations ( |Moe et al.|2009||Arav et al.|2012l ). [Booth fc Schaye 
( 2009 1 and Choi et al. ( 2012 \ have carried out careful studies 
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of the underlying numerical and physical assumptions that 
go into this feedback model, the latter paper finding that 
including the radial momentum of the AGN ejecta in the 
calculation dramatically increases the efficiency of mechan- 
ical feedback. 

Dubois et al. (f2010) developed a model of AGN involv- 



ing anisotropic injection of kinetic energy on scales simi- 
lar to the resolution of the simulation. They implemented 
this model in the adaptive-mesh-refinement code RAMSES 
( Teyssier|2002 1 and studied the implications of this feedback 
model on the properties of black holes and galaxy clusters in 
cosmological simulations a resolution of approximately one 
kiloparsec. 



Gaspari et al. (2012 I performed simulations of the evo- 



lution of the ISM of an elliptical galaxy in the presence of 
AGN feedback similar to that of Dubois et al. (20101. In 



terms of the galaxy model and the treatment of stellar evolu- 
tion and mass loss, the simulations are substantially similar 
to those presented here. The differences include the fact that 



the Gaspari et al. (20121 simulations are three dimensional 



rather than two dimensional, have much lower spatial resolu- 
tion (~ 150 pc rather than 0.25 pc, and include feedback only 
via anisotropic injection of kinetic energy, rather than the 
more comprehensive feedback model presented here. [Gaspari] 
et al. (20121 found that mechanical feedback efficiencies in 



the range of those observed by Moe et al. ( 2009 1 and Arav 
|et al.| ( |2012| | were sufficient to limit black hole growth, bal- 
ance atomic cooling, and lead to observationally reasonable 
gas density and temperature profiles, in substantial agree- 
ment with the conclusions of [Novak et al.| ( !2011j , ^Ciotti et al.| 
(20101, and other studies. 



The observational studies of Moe et al. ( 2009 1 and Arav 
|et al.| ( |2012| ) indicate that this form of mechanical feedback 
only uses a small fraction of the available energy. This mo- 
tivates the study of other forms of radiative feedback. 

2.2 Radiative momentum feedback 

It is well-known that the photons couple to the ISM weakly 
if only the small Thompson cross section is considered (e.g. 
Ciotti fc Ostriker|2001 1. Of the momentum available, only a 
fraction tes (the optical depth to electron scattering) cou- 
ples to the ISM, while there is no additional energy trans- 
fer to the gas in the case of pure Thompson scattering 
of photons having E <^ 1 MeV. Only a small fraction of 
AGN are known to be optically thick to Compton scat- 
tering, although studies of the X-ray background indicate 
that Compton-thick AGN may be common at high redshift 
( [Daddi et al.|200"7|. X-rays coupl e more strongly than lower 
energy photons ( Sazonov et al. 2005 1 due to both inverse 



Compton effects and coupling with resonant metal lines. 
Nevertheless, the Mbh cr relation ([Gebhardt et al.|2000[ 



Ferrarese & Merritt 20001 has been construed as evidence 



that black holes regulate their growth via momentum-driven 
winds ( |King|2003| ). 

Radiation also couples to the ISM through scatter- 
ing/absorption by dust grains as well as atomic resonance 
lines. The cross section for dust absorption in the UV for 
solar metallicity gas with a dust/gas ratio typical of the 
Milky Way can be 3500 times the electron scattering cross 
section ( jDrainepOOS l, while that for resonance lines can be 
the same or greater ( Proga et al.|2000 Sazonov et al.|2(305 1 . 



This means that radiative feedback due to line or dust ab- 
sorption can provide more energy and momentum than me- 
chanical feedback for column densities as low as 10~^** to 
10"^^gcm-^ 

The potential importance of radiative momentum feed- 
back has been known for some time (Haehnelt 1995: Ciotti fc| 
^Ostriker_2007| hereafter iCOTj), but the particular importance 
of dust in the context of AGN and starbursts was pointed 



out by Murray et al. (20051 and Thompson et al. (20051 



Simulations by Debuhr et al. (2010 2011 20121 model the 



effect of radiation pressure acting on dust at radii large com- 
pared to the black hole but small compared to the galaxy. 
Note that while the physical process mediating the interac- 
tion between the SMBH and the ISM is different from the 



mechanical model of Di Matteo et al. (20051, the photon 



momentum is converted to hydrodynamic motion on scales 
smaller than the resolution of the simulation (~ 100 pc). 

2.2.1 Simulations of radiative feedback 



Sazonov et al. ( 2005 1 investigated the important effects of 
X-ray heating and [Ciotti fc Ostriker| and collaborators have 
incorporated all of the standard electromagnetic processes 
into their one-dimensional treatments of the cooling flow ini- 
tiated AGN outb ursts ( |C07l [Ciotti et al.|2009[ [2010l [Ciotti[ 
fc Ostriker[[2012 K These investigations represent a different 
viewpoint from the "conventional" view that AGN are as- 
sociated with mergers. Here, the AGN activity attributed 
to the reprocessing of gas lost by the evolving stellar pop- 
ulation of an isolated and "passively" evolving early-type 
galaxy. In section [2.5[ we comment on recent observations 



that bear on this point. Recently Hensley et al. (20121 have 



added an improved treatment of dust creation and destruc- 
tion as well as a calculation of the dust temperature to these 
one-dimensional models. 

Feedback in the [Debuhr et al.| simulations is imple- 
mented via a force on the gas assuming that there is a fixed 
optical depth to IR photon scattering by dust at all times 
(r ~ 10). No self-consistent solution for the radiation field 
is sought. The radiative force is applied over one resolu- 
tion element near the black hole, so, although the physical 
model posits a radiative origin for the computed forces, in 
implementation their algorithm bears some resemblance to 
mechanical feedback schemes in that it is purely local and 
the photons do not exert forces over macroscopic distances 
in the simulation. 



Novak et al. (20111 used the same physical model de- 



scribed in Sazonov et al. ( 2005 1 and C07 (with the exception 



of the radiative forces on dust) to run two-dimensional sim- 
ulations of SMBH fueling and AGN feedback. The present 



work extends Novak et al. (2011 1 to include radiative forces 



on dust grains. We seek to carry out the radiative transfer 
calculation due to photon absorption and scattering by dust 
in a self-consistent fashion. 



Hambrick et al. (20111 and Kim et al. (20111 carried 



out simulations of black hole fueling and AGN feedback 
during galaxy mergers and found that the addition of X- 
ray heating by the AGN allowed the black hole to regulate 
its own growth by keeping its immediate vicinity hot with- 
out necessarily heating the majority of the gas in the galaxy. 
Hambrick et al. (2011 1 used cosmological smoothed-particle- 



hydrodynamics (SPH) simulations with a simple treatment 
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of X-ray photons from the AGN. Kim et al. (20111 per 



formed adaptive-mesh simulations of gala^cy mergers reach- 
ing a resolution of 15 pc. The simulations tracked mechanical 
feedback in the form of a coUimated outflow generated near 
the black hole as well as radiative feedback in the form of 
X-rays generated by the AGN. The X-rays transferred en- 
ergy and momentum to the surrounding gas, via by a three- 
dimensional Monte Carlo ray tracing algorithm. 



2.3 Angular momentum transport 

Several efforts have focused on understanding angular mo- 
mentum transport in the absence of feedback. This eases one 
of the severe computational constraints in the problem, since 
it is feedback that generates high gas temperatures and large 
gas velocities, necessitating very short simulation time steps 
and therefore very large computational costs per simulation. 



Levine et al. ( 2008 1 used adaptive-mesh cosmological simu- 



lations to examine gas transport from cosmological scales to 



1 pc from the black hole. Hopkins & Quataort (20101 car- 



ried out nested zoom-in simulations of the central regions of 
a galaxy undergoing a major merger. Both groups concluded 
that the mutual gravitational torques exerted by clumps of 
gas were sufficient to transfer angular momentum away from 
the central regions and permit accretion, and that the scale- 
free nature of gravity allowed the process to proceed in a 
nested fashion through all scales to the presumed black hole 
accretion disk. It is very attractive to split the problem of 
black hole fueling and AGN feedback into a fueling part and 
a feedback part. However, one would expect gas dumpiness 
to be dramatically reduced if AGN are able to heat the ISM 
significantly, which would in turn shut down the angular mo- 
mentum transfer process that these authors envision. If the 
gas clumps are shielded from the AGN photons by a cen- 
tral dusty torus, or if the clumps of gas are self-shielding, 
then the process may continue to operate in spite of AGN 
feedback. 



2.4 Sub-parsec simulations 

As high-resolution simulation efforts that resolve the Bondi 
radius become more common, technical details of the nu- 



merical solution become important. Barai et al. (2011 1 have 



performed a detailed study of the numerical and physical 
aspects of simulating Bondi accretion using SPH in the pres- 
ence of heating by a central black hole over scales of 0.1 to 
200 pc. They conclude that typical formulations of the artifi- 
cial viscosity term typically used in SPH can cause excessive 
heating near the inner boundary of the simulation. 



Dorodnitsyn et al. (20121 have studied the radiation 
hydrodynamics of the inner few parsecs of galaxies in the 
presence of dust in the the flux-limited diffusion approxi- 
mation. They show that the dusty torus around AGN is 
plausibly supported by radiation pressure on dust in an out- 
flowing wind. They are interested in smaller physical scales 
than presently concern us, and they are only interested in 
the limit of large optical depths to scattering by dust in the 
infrared. However, there is significant overlap in the basic 
physics between their work and the present paper. 



2.5 Observations concerning mergers, star 
formation, and AGN 

The observational situation regarding the links between 
black hole growth, star formation, and mergers has been 
somewhat confused, although it seems that a preponder- 
ance of evidence is now indicating that, although mergers 
may trigger AGN, the majority of AGN are not triggered by 
mergers. 



Pierce et al. ( 2007 1 used a combination of X-ray data. 



space-based optical data, and ground-based optical data, all 
part of the AEGIS survey ( [Davis et "ar]|2007[ ), to conclude 
that X-ray selected AGN preferentially reside in early-type 
galaxies and that although the X-ray selected AGN were 
more often members of close pairs, the companion was usu- 
ally undisturbed. This was interpreted as evidence against 
interactions triggering AGN. 

Using Sloan Digital Sky Survey data, Cistcrna s et al.| 
(2011| found a strong correlation between close galaxy pairs 
and galactic star formation, but no correlation between be- 



tween close pairs and AGN activity. Schawinski et al. ( 2010 1 



conducted a similar study and found superficially similar 
correlations, but based on their results connecting galaxy 



color to stellar population age (Schawinski et al. 20071 ar 



gued that the lack of AGN with close companions was due 
to a delay between the final merger and the onset of AGN 
activity. 



Ellison et al. (20111 have found that close pairs tend 



to have AGN activity, but that AGN do not have an el- 
evated close companion rate. They argue that mergers do 
indeed cause AGN activity, but that the majority of AGN 
are caused by other, presumably secular processes. 



Diamond-Stanic fc Rieke| ( |2012[ ) found that AGN activ- 
ity is correlated with nuclear star formation (within ~100 
pc), but not with global star formation on kiloparsec scales. 

^Rosario et al. ( 2011 1 used space-based optical and x-ray 



data taken as part of the CANDELS survey ( Grogin et al 



20111 to conclude that there is little difference between X- 



ray selected AGN and quiescent galaxies in terms of colors 
and stellar populations out to a redshift of three, arguing 
against a scenario where a significant fraction of AGN are 
triggered by mergers. As part of the same survey, IKocevski] 
et al. (20121 concluded that X-ray selected AGN are not 



morphologically different from a mass-matched sample of 
quiescent galaxies. 



3 THE RADIATIVE TRANSFER EQUATION 

In this section, we develop a method of treating the radiative 
forces forces on the ISM gas in a galaxy. Our goal is to arrive 
at an algorithm that gives the correct behavior in the many 
asymptotic limits required by the physical problem at hand. 

First we discuss the astrophysical requirements that 
motivate some of the choices we have made in defining the 
method. In the present work we are primarily concerned 
with radiative forces on dust grains, but the method works 
for any source of isotropic scattering or absorption opacity. 

More precisely, we derive three related algorithms for 
solving the radiative transfer equation. All of them arise 
from developing a set of moment equations for the specific 
intensity. For each method, we derive set of differential equa- 
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tions, the solution of which give the mean intensity and net 
flux of the photon field. 



The first method (Section 3.2 I gives a boundary-value 
problem where the physical content of the set of differen- 
tial equations is quite transparent. However, closing the set 
of differential equations in this method requires us to rely 
rather heavily on our physical intuition about the solution. 

The second method (Appendix|A| also gives boundary- 
value problem and is more aesthetically pleasing in a formal 
mathematical sense because the closure relation is derived 
along with the set of differential equations — it does not re- 
quire a physically motivated guess. The cost an additional 
differential equation to solve and the fact that the physical 
content of the equations is somewhat more opaque compared 
to the first method. 

The third method (Appendix p| is an attempt to cap- 
ture the essential behavior of the system while reformulating 
it as an initial value problem. The reason for this is that it 
is vastly simpler to efBciently and robustly obtain the nu- 
merical solution to an initial value problem. 

Note that in this work we use the first of these three 
methods and present the other two for reference and for 
future work. All three of the algorithms are self-consistent 
in the sense that we first solve for the photon field and then 
learn from this solution whether the radiation field is nearly 
isotropic (like the interior of a star) or highly directed (like a 
point source), whether the energy is carried by UV, Optical, 
or IR photons, and so on. We do not assume a priori that the 
system is in any particular asymptotic limit of the radiative 
transfer equation. 

We neglect angular force terms and although we use 
spherical symmetry to derive the differential equations, we 
argue in Section |3.2.f| that the method can be used with 
good results even in the case where the system is not spher- 
ically symmetric. 

3.1 Astrophysical preliminaries 

Over their lifetimes, galaxies explore essentially all of the 
asymptotic limits of the radiative transfer equation. Most 
of the time, most lines of sight are optically thin to radi- 
ation and the central SMBIf is not accreting significantly, 
so the radiation field within the galaxy is nearly isotropic 
and produced by a spatially distributed source (the stars 
themselves). However, most galaxies undergo brief periods 
of intense SMBH accretion, during which time the radiation 
field is dominated by the central point source and is not 
isotropic. During these times, most lines of sight remain op- 
tically thin, but a few lines of sight pass through the central 
dusty torus which is highly optically thick to UV/optical 
photons. Furthermore, many galaxies probably have spent 
time in a LIRG/ULIRG state involving intense bursts of 
star formation with a nearly spherical, highly optically thick 
ISM in the central regions. In this case, the radiation field 
may be sourced primarily by the central point source or 
the distributed stars, and in either case repeated scattering 
and absorption of UV and optical photons make the radia- 
tion field tend toward isotropy. Overall, the galaxy does not 
spend much time in this highly optically thick state, but in- 
teresting and important things are happening during that 
time. 

It is crucial to note here that all of the energy ab- 



sorbed in the UV/optical bands must be re-emitted as long- 
wavelength photons, because to neglect the re-emission is 
to ignore energy and momentum conservation. In particu- 
lar, ULIRGs are optically thick even to electron scattering 
over essentially all lines of sight, so the diffusion of the re- 
processed IR photons can have a significant impact on the 



state of the gas in the galaxy ( Murray et al. 2005 Thompson 
et al.|2005| ) 



In the optically thin case, the radiation field is highly di- 
rected in the region where the point source dominates. This 
region always exists at sufficiently small radii for a true point 
source. If and when the stellar light starts to dominate the 
radiation field, it becomes nearly isotropic for radii within 
the region that is producing the bulk of the photons. For 
significantly larger radii, the radiation field again becomes 
directed as the photons emerge from the galaxy on nearly 
radial paths heading toward infinity. 

If the ISM is optically thick to scattering, the radiation 
field quickly becomes isotropic even if the photons come 
from a central point source. In the UV and optical part 
of the spectrum, the opacities to scattering and absorption 
due to dust are of the same order of magnitude. Typically 
the scattering opacity is less than the absorption opacity 
(Draine'2003') , in which case the primary effect of scattering 
is to make the radiation field tend toward isotropy, rather 
than trapping photons so that they diffuse out of the galaxy. 
Therefore in adopting a closure relation we will assume that 
the radiation field becomes nearly isotropic (although still 
maintaining a net outward flux) when the optical depth to 
scattering or absorption is greater than unity. If the scatter- 
ing opacity is much smaller than the absorption opacity, this 
is incorrect: a highly directed radiation field from a point 
source is always highly directed even after many (absorp- 
tion) optical depths. However, for realistic dust properties 
the scattering and absorption opacities are of the same order 
of magnitude. 

If a UV or optical photon is absorbed rather than scat- 
tered, then the energy is reprocessed into the IR, but the 
energy is lost to the UV. Thus absorption will not tend to 
change the angular character of the radiation field (directed 
versus nearly isotropic). However, for declining density dis- 
tributions, the column density f pKdr is dominated by what 
is occurring at small radii, whereas the energy injected by 
stars, J e47rr^ dr (where e is the emissivity for stellar pho- 
tons) , is dominated by what happens at large radii unless the 
stellar distribution falls off as r~^ or faster. Therefore, even 
if absorption dominates over scattering (as it does in the 
UV), if the ISM is optically thick, the central point source 
will quickly be diminished compared to the radiation pro- 
duced by stars. 

Radiative transfer is difficult in general, in part because 
of the high dimensionality of the problem. The quantities of 
interest are functions of spatial position, photon propagation 
direction, frequency, and time, giving seven dimensions in to- 
tal. Furthermore, the solution is subject to non-local effects: 
photons tend to leave a system via optically thin "windows" 
if they are available. A nearly spherical gas cloud with op- 
tical depth r along most lines of sight and a few clear lines 
of sight comprising solid angle Q, will trap photons roughly 
as though it had effective optical depth Toff = min(r, il/47r). 
That is, in order to effectively trap photons within an opti- 
cally thick cloud, it is necessary to have a covering fraction 
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of order 1/r or else the photons wiU leak out of optically 
thin "holes." 



3.2 Mathematical treatment 

Using the physical insight provided by the discussion in 
Section |3.1| we now seek a mathematical formulation 
of the problem that is not overwhelmingly computation- 
ally intensive and is correct (or nearly so) in each of 
the identified limits (optically thin/thick, radiation field 
nearly isotropic/highly directed, dominated by distributed 
source/point source, and gas distribution nearly spheri- 
cal/highly non-spherical). 

We must distinguish the variables specifying the pho- 
ton propagation direction from spatial position; to this end 
we stipulate that 6 and (/> refer to spatial position while 6 
and 4> refer to photon propagation direction. Without loss 
of generality, in the case of isotropic scattering, as we con- 
sider in our treatment, the frequency-integrated equation of 
radiative transfer is 



dl{^., (t>) 
ds 



— —p{Ha + Hs)I + pK-sJ + e , 



(1) 



where s is the distance along the arbitrary direction of the 
ray, /i — cos 9, 9 is the angle between the ray propaga- 
tion direction and the z axis, <j) is the azimuthal angle, 
/ = dE/dt dfl dA is the specific intensity, J = J I djj. d(j>/4-K 
is the mean intensity, p is the gas density, Ka is the cross 
section for absorption per unit mass, Kg is the cross section 
for scattering per unit mass, and finally e is the emissivity 
per unit solid angle | Chandrasekhar 19601. All quantities 



are functions of spatial position and time, suppressed above 
for brevity. 

If we further adopt spherical symmetry, i.e., all quan- 
tities are independent of spatial position 9 and (p, then eq. 
(fTl) becomes 



dl 1 

t^TT ^ — 
or 



dp 



= —p{K,a + K.s)I + l^sJ + e . 



(2) 



Here we seek a solution to this equation by taking mo- 
ments in p and integrating away the dependence of the pho- 
ton field on propagation direction. The final ingredient is a 
physically reasonable closure relation to terminate the set of 
moment equations, to be discussed below. 

Note that although we write the radiative transfer equa- 
tion in spherical form, the numerical method we outline does 
not require spherical symmetry. The gas density, opacities, 
and photon field may vary as a function of 9 and cj). Further 
details including the method by which this is accomplished 
as well as advantages and disadvantages of our formulation 
are discussed in Section F3.2.1l 

The considerations in Section 13.11 can be summarized 
by the ansatz that the specific intensity to be composed of 
three terms: one isotropic, one mildly anisotropic, and one 
highly directed: 



/ = A + B/i + D5{^ - 1) , 



(3) 



where 5 is the Dirac delta function and the three coefficients 
(which are functions of time and radial coordinate) are ad- 
justed as described below. Incorporating the last term (with 
-D > 0) permits us to allow for phases when the central 



point-like quasar dominates the radiation field and to cor- 
rectly treat the photon field far outside of the galaxy where 
the specific intensity resembles that of a point source. Note 
that if we were to take D = in eq. ([3|, then we would 
arrive at the Eddington approximation, which is commonly 
adopted in stellar interiors. Taking eq. ([3| as a generaliza- 
tion of the Eddington approximation allows us to interpolate 
smoothly between the limits of a central point source and 
an optically thick situation. 

Taking moments in p of the specific intensity to remove 
the explicit dependence on the ray propagation direction 
gives: 



J: 



F : 



and 



1 

47r 



d(j> I I{fj,,(l>)dfi ^ A + 



d<t) 



D 



4ttB 
P Hp, 0) dp = — ^ -I- ttD , 



P^I{p^4>)dp = 



47rA ttD 
3c c 



(4) 



(5) 



(6) 



where J is the mean intensity, F is the energy fiux, and P 
is the radiation pressure. 

We must find some physically plausible way to decide 
whether the radiation field is nearly isotropic (in which case 
the standard Eddington closure relation P — it/3 applies, 
where P is the radiation pressure and u is the energy den- 
sity) or highly directed (in which case the flux is simply 
related to the zeroth and second moments of the specific 
intensity) . 

To find the appropriate equations, we take two moments 
of eq. (pl in two different cases. First, when the radiation 
field is mildly anisotropic {D — 0) with the closure rela- 
tion P = 47rJ/3c, yielding the differential equations for the 
photon field in the classic Eddington approximation: 

dL 



dr 

dj _ 3p{K,a + Hs)L 

dr 



(7) 



(8) 



167r^r^ 

Second, when the radiation field is highly directed {A — B = 
0,D — 1) with the closure relation P = 4nJ/c, appropriate 
for a point source. In this case the equation for dL/dr is the 
same and the equation for mean intensity becomes 

dJ _ 2J p(Ka+K-s)L 

dr r 



(9) 



IGvr^r^ 

The first term on the right-hand side gives the expected r~^ 
fall-off of the mean intensity if the radiation field is primarily 
directed (e.g. outside of the galaxy). 

We combine equations eqs. ([8| and ([9| for dj/dr by 
introducing the variable t that interpolates between the op- 
tically thin (radiation field like a point source) and optically 
thick (radiation field nearly isotropic) cases. This gives: 

dL . ,,^ . . ^^^^ 



= 4-Kr [E -4-KpKo.J) , 

2Ji (3 - 2t)p{Ka -f K,s)L 



(11) 



dr 

dJ 

dr r Wiv^r^ 

where t = when the radiation field is nearly isotropic and 
t = 1 when it is highly directed. In practice, we take i as a 
function of position to be: 
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l-f(r) 



where 



1 



p5(l- 



t) + e-"''(i-^in) 



Tin(r) = / p{Ka + Ka)dr, 

Jo 
and 



-Tout (r) 



p{K,a 



)dr. 



(12) 



(13) 



(14) 



Evaluating £ as r becomes very large or very small gives 
t = 1 at the inner boundary (near the AGN) and at the 
outer boundary (far from the galaxy). If the ISM is optically 
thick there will be some region in between where t is nearly 
0. This choice expresses our expectation about where the 
radiation field should be primarily directed versus primarily 



isotropic. The factors of five in eq. (121 are chosen to ensure 



that the value of t switches from the optically thin limit to 
the optically thick limit fairly quickly near optical depths 
of unity. The transition effectively occurs between optical 
depths of 0.8 and 1.2. 

Finally we specify the boundary conditions. As in the 
study of stellar interiors, our second order differential equa- 
tion for L{r) requires two boundary conditions, one at the 
center and one at the outer surface. Near the center of the 
simulation, the outgoing luminosity is only that provided by 
the central point source, the AGN. Far from the galaxy, the 
radiation field is expected to again look like a point source as 
scattering/absorption become negligible and all of the pho- 
tons leave the galaxy. The two natural boundary conditions 
are then 



-'^('"min) = LbH , 

and 



Jirn 



L{r 



A/ 



The force on the gas in a computational cell is 
_ L{r)dTdA 



(15) 



(16) 



(17) 



where dr is the optical depth (scattering plus absorption) of 
a given cell and dA is the area of the cell perpendicular to 
the radiation flux vector. 



3.2.1 Comments on the spherical symmetry assumption 

We are primarily interested in what happens to the ISM as 
a function of radius rather than angle. Angular effects are 
important and it is of crucial importance to break spher- 
ical symmetry so that hydrodynamic instabilities such as 
convection and the Rayleigh- Taylor instability can operate. 
However, we are concerned mainly with whether gas makes 
it to small radii (and accretes onto the black hole) or to large 
radius (and escapes the galaxy), not where on the sphere the 
gas enters the SMBH or leave the galaxy. Thus we neglect all 
angular radiative forces and formulate the radiative transfer 
as occurring exclusively along rays from the center of the 
galaxy. 

Similar approximations have been used in numerical 
simulations of photon and neutrino transport in super- 
novae for some time (e.g. Brandt et al.|2011 and references 



therein). One approach is to spherically average the matter 
distribution, solve for the radiation field, and then calculate 



forces on matter via eq. (171. Another approach is to treat 



each ray as a separate problem, in which case all intensive 
quantities (gas density, emissivity of stellar photons, etc) are 
treated as though they obtained over the whole sphere. 

We found that the first approach has dramatic fail- 
ure modes where it wildly overestimates the outward radial 
forces on the gas in certain situations. As an example, con- 
sider an optically thick cloud covering a small solid angle 
in an otherwise optically thin medium. In reality, UV pho- 
tons deposit momentum in a thin layer at the surface of the 
cloud, and the total momentum imparted to the cloud de- 
pends on the solid angle of the cloud. However, spherically 
averaging the matter distribution may result in the ISM be- 
coming optically thin at the cloud radius, and the solution 
for the photon field will indicate that the UV flux is high 
throughout the cloud. In this case, essentially all of the dust 
grains in the cloud are able to absorb UV photons and the 
momentum imparted to the cloud will be wildly overesti- 
mated. 

Instead, treating each ray as a completely independent 
problem has much better behavior in the sense that when 
the algorithm makes an error in computing the force on a 
particular gas cell, fractional error tends to be of order unity 
rather than orders of magnitude. In the previous example, 
this method gives the correct force exerted by UV photons. 
For the IR photons, the method assumes that they must tra- 
verse the cloud on their way out of the galaxy, while in reality 
they would leave the system via the optically thin ambient 
medium. Therefore the method makes a a fractional error of 
order the IR optical depth when computing the force. 

One drawback of the latter approach is that if a par- 
ticular ray (say, near the pole) has deflcient star formation 
compared to the rest of the galaxy, then the forces exerted 
by the stellar photons will be anemic compared to the true 
solution, where stars from near the midplane emit photons 
that help push on gas near the poles. Therefore, we average 
source terms over a sphere before flnding the solution to the 
radiative transfer equations on each ray. 

Note that this method allows the radiation field to dif- 
fer in an unlimited way from one ray to the next. This is 
expected for example in the UV where the scattering opac- 
ity is less than the absorption opacity, and the interior of an 
individual "cloud" along a line of sight can be shielded from 
the AGN UV luminosity. However, in the IR we consider 
the absorption opacity to be zero, and it no longer makes 
sense that the radiation field could differ substantially on 
adjacent rays. Therefore we also spherically average the IR 
source terms (in particular IR generated by absorption of 
UV or optical) before solving for the luminosity as a func- 
tion of radius in the IR. 



3.3 Numerical implementation 

The numerical solution of the radiative transfer equation is 
complicated by the fact that one boundary condition is spec- 
ified at r-min while the other is specified at Tmax (note eqs. 



(15 1 and (16l). Thus, integration from either side requires 



a guess that must be refined. This is identical to the situ- 
ation in treating the radiative transfer in stellar interiors. 
Conceptually, this is simply because we know the central lu- 
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minosity, but we do know know the intensity of the radiation 
field near the center. If the scattering opacity is large, pho- 
tons will build up in the galaxy forming a nearly isotropic 
radiation field until the gradient of the photon density is 
sufficient to drive photon diffusion out of the galaxy at the 
rate required to carry the specified luminosity. 

In fact, the numerical solution to the radiative transfer 
equation in stellar interiors is considerably simpler than in 
the present problem because within stars, one can exploit 
the fact that the plasma is everywhere close to local ther- 
modynamic equilibrium. The plasma is assumed to quickly 
thermalize any energy sources or sinks so that the plasma 
is well-described by a single temperature at any point. In 
galaxies, the photon field is not close the thermodynamic 
equilibrium and we must track energy as it is transferred 
from x-rays to UV photons and eventually to IR photons. 
Treating stellar interiors in an analogous way would mean 
tracking each energy band of photons separately as the star's 
luminosity was transported through the plasma, from the 
primary MeV photons produced by fusion reactions, to the 
UV photons that eventually emerge from the photosphere. 

To find the solution of eqs. \1Q\ a nd ( |11[) s ubject to the 
correct boundary conditions eqs. (151 and (16 1, the "shoot- 
ing method" with Newton- Raphson iteration (Press et al 



|1992[ ) works well up to optical depths of a few. Beyond that, 
the absorbing gas acts as a very effective "screen" that hides 
conditions near the center of the simulation from those at 
the outer edge. This makes it very difficult to find the cor- 
rect solution matching both boundary conditions at once. 
We would like to simulate AGN feedback in galaxies with 
gas surface density up to Compton-thick, implying an opti- 
cal depth to dust absorption in the UV of several thousand. 
The shooting method is too unstable to use in this case. 

It is true that given such a large optical depth, very 
little of the primary UV radiation escapes the galaxy; it 
comes out in the IR. However, it is incorrect to discard the 
primary UV radiation completely, or to immediately transfer 
it to the IR. While this would correctly give the emergent 
energy, the UV radiation can have a very large effect near 
the center of the galaxy before it is transferred to the IR. 
For for solar composition and dust/gas ratios, they can limit 
the accretion rate to 3 x 10~'*LEdd- 

At large optical depths, the more powerful "relaxation" 
iterative method ( [Press et al.|1992| ) is required to adjust the 
solution over the entire domain of integration while ensuring 
that the proper boundary conditions are met at both ends of 
the interval. Newton's method is used to iterate over trial so- 
lutions to arrive at the solution to the differential equations, 
where the independent variables of the Newton iteration are 
taken to be the values of L and J at each spatial point. 
The method takes the coupling between variables at neigh- 
boring points into account, allowing the entire solution to 
be updated without overcompensating (which the shooting 
method is prone to do). Given a good initial guess for the 
iteration, the relaxation method finds a reasonably accurate 
solution even for very large optical depths. The initial guess 
is conveniently provided by the solution that pertained dur- 
ing the previous simulation time step. 

Two additional numerical schemes are defined in appen- 
dicies. The method described in Appendix [A] allows one to 
avoid using the variable t and its somewhat ad hoc definition 



equation to solve and a somewhat more opaque relationship 
between the equations and their physical content. The algo- 
rithm described in AppendixlBlis an attempt to preserve the 
basic physics of the method while transforming the system of 
differential equations from a boundary value problem to an 
initial value problem. This allows one to use a substantially 
simpler numerical method to obtain the solution. 



4 HYDRODYNAMICAL SIMULATIONS 

For the full simulations, we will use the 2D numerical code 
described in Novak et al. (20111. The code implements 



sources and sinks of mass, energy and momentum relevant 
for an early-type galaxy including stellar mass loss from 
planetary nebulae. Type la and Type II supernovae. A full 
description of the input physics can be found in |Novak et al.| 
( |2011[ ) and |Giotti fc Ostriker] ( |2012[ ). We describe the sim- 
ple treatment of dust creation and destruction used in the 



present work. We refer readers to Hensley et al. (20121 for 
a more complete implementation of dust physics applied to 
one dimensional models. 



4.1 Timescales 

The light crossing time of the galaxy is by far the shortest 
timescale in the system and therefore if the ISM is optically 
thin, photons are considered to act instantaneously for the 
purpose of our simulation. When the ISM becomes optically 
thick, the radiative diffusion time is 



^diff 



Rt 



(18) 



where r is the optical depth. Typically the temperature of 
the ISM is not far from the virial temperature of the galaxy, 
so that the dynamical time and the sound crossing time are 
nearly equal. The latter is 



ts 



R 



(19) 



Thus the radiative diffusion time will be shorter than the 
dynamical time and sound crossing time as long as r Sj 
c/cs ~ 1000. For our purpose radiative diffusion takes place 
primarily in the IR, and the maximum optical depth in the 
IR is of order 10. Therefore the radiation field is considered 
to reach equilibrium instantaneously, and radiative forces 
are considered to act instantaneously. That is to say that 
we solve all of the radiation transfer equations assuming a 
steady state. 



4.2 Photon energy bands 

We divide the radiative energy output of the AGN into broad 
bands where the relevant physics is significantly different: IR 
(below 1 eV), Optical (for our purposes, non-ionizing: 1-13 
eV), UV including soft X-rays where the absorption cross- 
section is greater than the Thompson cross-section (13-2000 
eV), and hard X-rays where the interactions are dominated 
by Compton scattering (^2 keV). The assumed division of 
the AGN output into the different bands is given in Table 



m 



in eq. ( 12 1 at a cost of one additional additional differential 



The radiative transport in the X-ray band is relatively 
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Table 1. Assumed unabsorbed QSO energy output by band 
( [Sazonov et al.|2004^ and dust opacities for solar metallicity gas 
with Milky Way dust-to-gas ratio expressed in terms of the elec- 
tron scattering opacity, kes = 0.4 cm-^ g~^ 



tions given by eqs. (10 1 and (111. The boundary conditions 



Band Energy 



Fraction Dust Opacity 



X-ray E >2 keV 0.1 

UV 13 eV <E <2 keV 0.35 3430 kes 

Optical 1 <E <13 eV 0.25 857 kes 

IR E <1 eV 0.3 5.71 kes 



simple. Few AGN are observed to be Compton thick and 
our simulations rarely reach such a high column density, so 
the ISM is always treated as optically thin. X-rays dominate 
the heating of the gas to temperatures of order the galac- 
tic virial temperature. We treat heating and cooling due 
to photoionization and Compton heating/cooling by atomic 
lines via the formulae given in Sazonov et al. (20051. X- 



rays also affect the gas momentum via Compton scattering 
and photoionization. Compton scattering by X-rays is han- 
dled along with the bolometric AGN luminosity owing to the 
wavelength independence of the Compton cross section. Mo- 
mentum transferred to the ISM by photoionization is taken 
to be 



dp _ H 

dt c 



(20) 



where H is the photoionization heating rate. It should be 
noted that this is a lower limit to the total effect of pho- 
toionization on the momentum of the ISM. The Sazonov 
formulae give the photo heating rate due to X-rays, but ne- 
glect the heating at lower energies due, e.g. to optical and 
UV photons. 

An exact treatment of the effect of optical and UV pho- 
tons necessitates tracking the temperature as well as the ion- 
ization state of the gas — this quickly becomes prohibitively 
complex. Instead we note that X-rays are required to raise 
the gas temperature to the same order as the virial temper- 
ature for a massive galaxy, and therefore drive an outflow 
into the IGM. Lower energy photons may affect the detailed 
state of the gas in the galaxy, but they cannot affect whether 
or not the gas remains bound to the galaxy. Compton heat- 
ing can be important if the Compton temperature of a given 
object is large. However, for typical objects Compton heat- 
ing and cooling are sub-dominant. While optical and UV 
photons acting on atomic resonance lines are unable to raise 
the gas temperature to the galactic virial temperature, they 
nevertheless contribute to the force on the ISM. For the time 
being we neglect them. 

Radiative transport in the IR is similarly simple. At 
long wavelengths, the dust absorption opacity is much 
greater than the scattering opacity ( Draine|2003 |. However, 
the absorbed energy is re-radiated by the dust grains and 
thus the process is more similar to scattering than true ab- 
sorption. In order to conserve energy, we simply treat the 
process as scattering. Since no energy is ever lost from the 
IR band and must leave the galaxy eventually, the IR lu- 
minosity as a function of radius is simply the integral of all 
source terms within the radius. 

The UV and optical bands are considerably more com- 
plex. Each band satisfies a separate set of differential equa- 



and method of solution are as described in Sections 13.21 and 
|3.3| applied to each band separately. Radiation absorbed in 
both bands contributes to the IR. 



4.3 Opacities 

At a minimum, photons emitted from the stars and AGN 
impart momentum to the ISM via Compton scattering. The 
Compton cross section is wavelength independent and coher- 
ent up to energies of order 511 keV, making the effect easy 
to treat numerically. The force on a cell is given by eq. ( 17l 
where L is the bolometric luminosity of stars and the AGN 
within the radius r. This is correct even when the ISM is 
Compton thick because the process involves only scattering 
and no absorption (up to energies approaching 511 keV). 

The effect of absorption by dust is more complicated 
because the opacities depend strongly on wavelength, and 
hence it is necessary to track the luminosity as a function 
of photon energy and radius. The dust opacities for each 
band for solar metallicity metallicity gas with a Milky Way 
dust-to-gas ratio are given in Table [l] 

Finally, outgoing photons can be absorbed by resonant 
atomic lines. In the case that the atom enters an excited 
state, the photon produced by the decay of the excited state 
is emitted isotropically. Thus the photons impart momen- 
tum, but no net energy to the gas. These bound-bound 
transitions can have very large opacities, but only over a 
narrow range of frequencies, so we neglect this effect in the 
present work. In the case of photoionization, the energy in 
excess of the binding energy goes into heating the gas, im- 
parting both energy and momentum to the gas. For optical 
and UV photons, accurate calculation of the photoionization 
opacity requires detailed knowledge of the temperature and 
ionization state of ISM gas that the present simulations lack. 
Therefore we leave the inclusion of photoionization opacity 
in the optical and UV bands for future work. For the high- 
energy photons that dominate the heating of the gas (see the 
discussion in the previous section and Sazonov et al.|2004 |, 
the photoionization force is included through eq. (|20 
discussed above. 



4.4 Dust creation and destruction 

Dust in ETGs is thought to be produced primarily in plane- 
tary nebulae and asymptotic giant branch outflows, and de- 



stroyed in hot gas by sputtering. Draine & Salpeter ( 1979 1 



give a theoretical estimate of the dust destruction time as a 
function of grain size: 



U. 



cstroy 



10' 



yrao.m 



I'a+T,-')., 



(21) 



where ao.i — a/0.1 fim, ni = p/mp cm~^, and Te = 
T/10^ K. Dust is also commonly observed to be produced 
via supernovae. Therefore, when gas is added to the com- 
putational grid by explicit stellar source terms, we also add 
dust in proportion. 

Draine ( 2009 1 has recently argued that dust grain 



growth is predominantly by gas-phase growth in the ISM 
rather than in supernovae. If this is the case, we should add 
an explicit dust source term that causes the dust to gas ratio 
to relax to a specifled value when the gas is cold and star 
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forming, regardless of whether all of the gas in the cell is 
actually processed through supernovae. 

In this scenario, dust grains are thought to be built up 
by accretion of single atoms from the gas phase, the inverse 
of the process that destroys grains at high temperature. The 
only difference is that grain growth requires that a dust grain 
collide with a metal atom rather than a more abundant hy- 
drogen atom. Metal atoms are both less abundant and move 
more slowly in thermal equilibrium. Therefore the timescale 
for creation is easily obtained by scaling from the destruc- 
tion timescale: 



Table 2. Parameters governing stellar luminosities (see text for 
definitions) 



Band Fraction 

euv 8.65 X 10-5 

eopt 1-24 X 10-3 

tuv 2.57 Myr 

tOpt 154 Myr 



= 10' yr A^^^Z-'ao.m^'T^^^il -f 10 Te) , 



(22) 



where A is the atomic weight of the atoms in question (here 
taken to be A = 12 for Carbon), Z = Pmotai/p is the metal- 
licity of the gas, and the (1 -I- 10 Te) term is added to ensure 
that dust creation takes place only at temperatures below 
lO' K in the simulation. 

The dust density in a given cell is represented by a 
passively advected tracer field with appropriate source terms 
for dust creation and destruction. The source and sink terms 
for dust in a given cell are 



Pdust 



Pdust pZ - 



' Pdust 



idc 



t'CTi 



(23) 



Dust grain growth terminates when all of the metals are in 
dust grains. In this paper the treatment of dust is signifi- 
cantly improved with respect to the basic approach in the 
series of papers by |Ciotti &: Ostriker] however, in |Hensley| 
|et al.| ( |2012[ ) , the dust treatment is discussed in full general- 
ity, albeit in the context of ID simulations. 

Including advection, the partial differential equation for 
the dust density is 



Qpdust 9(t^r pdust) 1 djVePduBt) 

dt dr r 89 



td. 



Pdust pZ — Pdust ,„.s 



estroy 



tc 



The radiative force on a cell due to absorption by dust 



is obtained from eq. ( 171 



AF- 



pD L{r)dTdA 
p Anr^c 



(25) 



where dr is the optical depth if there were no dust depletion, 
and the effect of dust depletion is taken into account by the 
density ratio. 

One possible improvement to our treatment of dust cre- 
ation and destruction would be to track the dust size distri- 
bution directly via moments. Our current implementation 
essentially assumes that all dust grains are maintained with 
a given fiducial size of 0.1 ^m and are then eroded via sput- 
tering or grown via metal atom deposition. 

4.5 Stellar source terms 

Since stars are sources of UV and optical light that con- 
tribute to the solution to the radiative transfer equation 
discussed here, we briefly review our implementation of star 
formation and generation of star light. Gas that forms stars 
also contributes to the energy balance via supernovae, our 
implementation of which is discussed in detail in |C07[ 

Gas is removed from the galaxy by star formation at 
the instantaneous rate given by the standard formula based 



on the Kennicutt ( 1998 1 law: 



Pgas — 



VtovmP 



Viovm = 0.1, 



3fcflT 



f-cool — 



2nA(r) ' 



Cdyn — IHini^t^jcansi i-rot J ; 



32Gp 



Crot — 



2Tvr 
Vc{r) ' 



(26) 

(27) 
(28) 

(29) 



where n^A(r) is the volumetric cooling function and Vc{r) 
is the circular velocity at radius r in the galajcy (Ciotti & 



Ostriker] ([20121). 

Mass that is transferred from gas to stars at radius r is 
first transferred to a "forming" state by convolving the star 
formation rate with an exponential kernel with a time con- 
stant of 3 Myr, designed to account for the finite minimum 
star formation given by the Kelvin-Helmholtz timescale as 
well as the minimum stellar lifetime given by the combi- 
nation of the Eddington luminosity and the nuclear energy 
available to a massive star: 



Pfo 



.w 



1 

iKH 



Pgas 



it') 



*KH dt 



(30) 



These stars also contribute photons to the UV/optical 
luminosity of the galaxy via convolution with a similar ex- 
ponential kernel: 



dLs 
~dV 



it) 



1 
t^ 



<^bc pgas(i')e '^ dt' 



(31) 



where B refers to band and takes the values UV or optical. 
Values for es and ts for each band are given in Table [2] 



4.6 Radiative cooling and heating 

Photons generated by coUisional excitation of atomic lines 
or free-free radiation can only cool the ISM insofar as the 
photons can actually escape. If the galaxy becomes optically 
thick to IR radiation, then these photons can contribute 
significantly to the forces experienced by gas in the ISM as 
they diffuse out of the galaxy. Therefore we add the cooling 
radiation to the distributed UV (optical) source terms if the 
local gas temperature is above (below) 10''' K. The radiative 
formulae are discussed in Ciotti & Ostriker (20121. 



4.7 Broad line wind 

Our simulations also include a treatment of the broad-line 
wind generated near the SMBH. This is implemented as an 
inflow term on the edge of the computational grid (Novak| 



et al. 2011 Ciotti & Ostriker 2012 1. Although the present 



paper focuses on radiative AGN feedback, we summarize 
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Figure 1. A snapshot of the A2 simulation during a quiescent phase. Upper left, log gas density in protons per cubic centimeter. Upper 
right, log sound speed. Lower left, radial velocity. Lower right, the log of the dust to gas ratio normalized to the value in the Milky Way. 
Gas that enters the ISM via the broad-linc wind from the AGN is assumed to be dust rich. At small radii, densities and temperatures 
are large enough to efficiently destroy the dust grains, resulting in dust suppression by a factor of up to 1000. 



our implementation of mechanical feedback because it has 
a large effect on the regulation of black hole growth. See 



Novak et al. (20111 and Ostriker et al. (20101 for further 



discussion of mechanical feedback. 

The BAL originating near the SMBH provides energy, 
momentum, and mass from a wind to the ISM according to 
the equations (Ostriker et al.|2010|: 



Mbh 



Mi, 



1 + 7? ' 



r\ = 



2ewc^ 



(32) 



Ew = em/A/bhc , 
and 



Pw = 



2ewc^MBu 

Vw 



Mw = 



2ewc^MBH 



(33) 



(34) 



where vw is the velocity of the BAL, taken here to be 10,000 
km s~^. 

Independently of the mechanical feedback model, the 
radiative luminosity of the AGN is given by 
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Figure 2. A snapshot of the A2 simulation during a black hole outburst. All gas near the center of the galaxy is depleted of dust by a 
factor of at least 30 and up to 1000. 



L — ^emMbhC 



(35) 



where the electromagnetic efficiency is given by the advec- 
tion dominated accretion flow inspired ( Narayan fc Yi|1994 1 
formula: 



eem = 



eoArh 
l + Arh' 



(36) 



and A = 100 and eo = 0.1. The dimensionless mass accretion 
rate is 



Mbh eoMsHC^ 



MEdd 



L. 



(37) 



where LEdd is the Eddington luminosity. 



Although our other work has considered models where 
ew is a function of the accretion rate, in the present work we 
assume that ew is a constant with values given in TablelS] In 
previous papers these are called A models, and we maintain 
that nomenclature here. Our AG model has the same value 



of ew as that adopted by Di Matteo et al. (20051, and our 
A2 model has a lower value, favored by recent observations 
( Moe et aL||2009[ |Arav et al.|2012 l. 
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Table 3. Values for the mechanical wind efSciency considered in 
this work. 



Model 



<^w 



AO 
A2 



5 X 10"^ 
10-* 



4.8 Galaxy model 



The total gravitational potential of the model galaxy is as- 
sumed to be a singular isothermal sphere plus a point mass 
for the central BH. This is good agreement with observations 



of the total mass profile of early-type galaxies ( Gavazzi et al 



2007 2008 1. For simplicity we maintain this model in to the 



smallest radii, although more complicated models may be 
more appropriate inside of a fraction of the half-light ra- 
dius. The velocity dispersion parameter of the isothermal 
potential is 260 km s~^. The gas is not self-gravitating. The 
stellar distribution is given by a Jaffe profile with a total 
mass of 3 X lO^^M© and a projected half-mass radius of 6.9 
kpc. The mass-to-light ratio is assumed to be spatially con- 
stant and is equal to 5.8 in solar units in the B band at 
the present time. The structural and dynamical properties 
of these galaxy models are discussed in detail in |Ciotti et al.| 
(20091). 



We note that although the galaxy model is spherically 
symmetric and the physics as implemented is either sper- 
ically symmetric or at least symmetric under inversions of 
the z axis, numerical noise is quickly amplified chaotically 
in the simulations so that the physical state does not remain 
up/down symmetric for very long. 



5 RESULTS 

All of the simulations presented here use the method de- 
scribed in Sections |3.2| through |3.3| to solve for the photon 
field, where the true boundary value problem is solved via 
the relaxation method. 

Figures [l] and |2] show two snapshots from the "fiducial" 
A2 simulation with mechanical AGN feedback efficiency of 
ew = 10~* and all of the physics thus far described enabled. 
The dust to gas ratio is initialized to the Milky Way value 
throughout the galaxy. This simulation uses the "standard" 
recipe for dust grain creation where dust is created in su- 
pernova and dust grains do not not grow via collisions with 
gas-phase metal atoms in the ISM. However, all gas that 
enters the simulation grid does so with a Milky Way dust- 
to-gas ratio. 

Figure[l]shows that near the center of the galaxy, dust is 
efficiently destroyed so that forces on gas due to absorption 
of photons by dust grains are not very effective at prevent- 
ing gas from falling into the center of the galaxy. At large 
radius, the temperatures are sufficient for sputtering to be 
effective, but the gas densities are low enough that the grain 
destruction timescale is long. 

In this picture, cold gas is generated by unstable cooling 
of the ISM. Cold gas loses pressure support and falls toward 
the center of the galaxy on a free-fall timescale. Cold gas 
is generated from hot gas, and the dust grains in the hot 
gas have been eroded by sputtering for for at least a cooling 



time. Therefore radiative forces on dust grains do not play 
a major role in the dynamics. Even if dust grains are re- 
plenished in the ISM without being processed through stars 



as suggested by Draine ( 2009 1 , the grain growth timescale 



is at best comparable to the infall timescale. Given their 
initially depleted state, dust grains must have many growth 
timescales available in order to replenish the dust. 

If radiative forces on dust are to play a role in black 
hole accretion, then the gas feeding the black hole must 
either never have been hot, or it must be cold for a long 
time before approaching the black hole. If the gas had been 
hot, it cannot return to the cold phase until a cooling time 
has passed, sufficient to significantly deplete the dust grains. 
Cold gas that is not rotationally supported falls to the cen- 
ter of the galaxy too quickly to build up significant dust. 
However, dust lanes are observed in a significant fraction 



of elliptical galaxies ( Hawarden et al. 111981 van Dokkum & 



Franx 1995 1. If the black hole is fed by gas processed through 



a rotationally supported galactic-scale disk, then dust grains 
would have a chance to build up while the gas spends a long 
time orbiting the black hole. 

Figures |3] and |4] show black hole growth over 1 Gyr for 
the AO and A2 simulations with mechanical feedback effi- 
ciencies of 5 X 10"'^ and 10"'* respectively. In the former 
case, dust has little effect on the dynamics of the simulation 
because the mechanical feedback dominates over everything 
else. In the latter case, dust makes some difference, although 
it is not a dramatic effect. Furthermore, Figure |4] shows 
that occasionally the simulations with different assumptions 
about dust arrive at very similar black hole masses in spite of 
the black hole masses having been different at earlier times 
(e.g. at 200 Myr). This indicates that sometimes the dust af- 
fects when gas reaches the black hole, but not i/gas reaches 
the black hole. 

In our simulations, dust is not negligible in determining 
black hole growth, but neither is it dominant. Extreme as- 
sumptions about the dust-to-gas ratio ranging from no dust 
at all to constant Milky Way ratios makes at best a factor 
of two difference in the mass added to the black hole. 

Figures [5] and |6] show the cumulative distribution of the 
ratio of the black hole accretion rate to the Eddington rate. 
The differences are not huge, but in the A2 case, dust opac- 
ity apparently does prevent the black hole from reaching the 
Eddington rate: the maximum is about 20% of the Edding- 
ton rate. The fraction of the time when the SMBH would 
be identified as a quasar [L/L-Edd > 0.1) is moderate. 

Figures [7] and Is] show the effect of various physical pro- 
cesses on the momentum balance of the ISM. These figures 
immediately show which physical processes dominate and 
which are negligible. It is important to note that these fig- 
ures show the time averaged force per unit volume, so the 
dominant forces could change in the midst intense bursts of 
AGN accretion and star formation. 

Most of the forces act on each grid cell separately, mak- 
ing it easy to tabulate how much momentum each separate 
process imparts to the ISM. Mechanical momentum is differ- 
ent: it is injected at the first grid cell and then the conserva- 
tion equations of hydrodynamics transport that momentum 
outward until it mixes with the ISM and the wind stops. 
However, the wind termination is not abrupt and there are 
significant radial motions in the ISM even in the absence of 
a mechanical wind from the AGN, so determining the exact 
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Figure 7. Time averaged radial momentum source and sink terms over 1 Gyr for the A2 simulation. Each line shows the time averaged 
force per unit volume imparted by each physical process as a function of radius. Curves have been multiplied by (r/kpc)^ to flatten them 
out. The thin rising blue line is gravity. Heavy lines are forces directed inward. Where a given physical process has both a heavy and 
light line, the force is sometimes directed inward and sometimes directed outward. The red line is for mechanical momentum. Mechanical 
momentum is injected at the first grid cell, so the quantity plotted is the force per unit volume if the momentum is deposited uniformly 
within a sphere of radius r. The green lines are for electron scattering of photons produced by the AGN or by stars, the yellow line is 
for photoionization by the AGN, the pink lines are for dust opacity for photons produced by the stars in different bands, and the blue 
lines are for dust opacity for photons produced by the AGN in different bands. The mechanical momentum accounts for the bulk of the 
momentum injected by feedback, and it will dominate at the radius at which the momentum is actually deposited into the ISM (typically 
~1 kpc). Photoionization by the AGN (AGN lines) is important inside a few hundred parsecs, outside of which AGN UV photons acting 
on dust are important. 



radius at which the mechanical wind deposits momentum is 
somewhat ambiguous. The quantity plotted as the red line 
in Figures [7] and [8] is the force per unit volume if the mo- 
mentum is deposited uniformly within a sphere of radius r. 
Since the lines are also multiplied by (r/kpc)^, mechanical 
momentum shows up as a straight horizontal line. At very 
small radii, the "dentist drill" effect ensures that not much 
radial momentum is actually mixed into the ISM. The me- 
chanical wind typically does not extend to the outer edge of 
the simulation grid, so likewise little momentum is actually 
deposited at large radius. Therefore the curve showing where 
the mechanical momentum is actually deposited would be 
low at small radii, low at large radii, nearly equal to the red 
line at intermediate radii (typically about 1 kpc) where the 



radial momentum is from the wind is actually mixed into 
the ISM. 

For the A2 simulation with dust destruction allowed, 
the most important sources of radial momentum are the 
mechanical wind, photoionization by AGN photons, and ab- 
sorption of UV AGN photons by dust grains. When dust de- 
struction is not allowed and therefore the effects of dust are 
maximized. Figure [S] shows that the most important sources 
of momentum are the mechanical wind, photoionization by 
AGN photons, and dust scattering and absorption of optical 
and IR photons originating from stars. In the latter, max- 
imal dust case, the effect of stellar photons on dust grains 
comes becomes more important than the effect of AGN pho- 
tons on dust grains for three reasons: 1) the black hole mass 
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Figure 8. Radial momentum source and sink terms over 1 Gyr for the A2 simulation, assuming a Milky Way dust-to-gas ratio at all 
times. The lines and labels are the same as Figure [7] except for the cyan line, which shows the very small amount of momentum removed 
by limiting the gas velocity to 30, 000 km s~^. The most important momentum sources are the mechanical wind, photoionization by the 
AGN, and scattering of stellar IR photons by dust grains within about 1 kpc. Outside of 1 kpc, absorption of optical stellar photons by 
dust grains is the dominant source of momentum. 



is about a factor of two smaller, so the black hole produces 
fewer photons; 2) about 20% more stars are formed, so stel- 
lar luminosity is increased; 3) the stars that form are much 
more centrally concentrated (about a factor of two more 
stars inside of 1 kpc), so the outgoing stellar photons see a 
larger optical depth to dust absorption because the density 
of the ISM falls steeply with radius 

It is somewhat surprising that dust seems to have such 
a small effect on the black hole growth in the simulation, 
given that the opacity in the UV due to dust is 3500 times 
the electron scattering opacity. Therefore, at least in the 
"Max Dust" case where a Milky Way dust to gas ratio is as- 
sumed for all gas, one might expect that the effective limit 
to the black hole accretion rate is ~l/3500 times the Ed- 
dington rate. However, the crucial crucial difference is that 
radiative transfer is governed by the absorption opacity in 
the UV while electron scattering and the transfer of thermal 
IR photons is a scattering problem. Once the UV has been 
absorbed, the radiative transfer in the IR is, again, essen- 
tially a scattering problem with kir = 5.71kes; however. 



dust to gas ratios in elliptical galaxies are depleted as com- 
pared to the Milky Way by more than the modest factor of 
five required for electron scattering to dominate over dust 
scattering. 



As noted by Thompson at al. ( 2005 I and Murray et al 



( 2005 1 in the context of this problem, the force exerted on 



a shell of gas by photons where scattering dominates is: 



/ 



c 



(38) 



where L and r are the luminosity of the source and optical 
depth of the shell, both defined in the given frequency band. 
The energy in photons cannot be destroyed, so if the shell 
is optically thick, then photons build up until the gradient 
of the photon number density is sufficient for diffusion to 
carry the input luminosity. This is why the Eddington lu- 
minosity effectively limits the black hole accretion rate even 
for Compton-thick sources, with the caveat that some en- 
ergy is transferred to thermal energy in the gas by inelastic 
scattering of photons with energies of order 1 MeV. The 
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Figure 3. Black hole growth over 1 Gyr in the AO simulation 
for different assumptions about the dust. The blue line shows the 
result of tracking dust density with a tracer field and accounting 
for dust creation and destruction. "No Dust" has all dust opacities 
set to zero. "Max Dust" assumes that all gas in the simulation has 
a Milky Way dust-to-gas ratio at all times. The high mechanical 
feedback efficiency in this simulation dominates over all other 
concerns and the presence or absence of dust does not make a 
significant difference. 
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Figure 5. Cumulative time (in units of total elapsed time) where 
the ratio of the black hole accretion rate to the Eddington rate 
is above the given value for the AO simulation with different as- 
sumptions about the dust to gas ratio. Dust does not make a big 
difference, although the "Max Dust" simulation spends somewhat 



more time with an Eddington ratio of 10 
the other simulations. 
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Figure 4. Black hole growth over 1 Gyr for the A2 simulation. 
The presence or absence of dust makes up to a factor of two 
difference in the black hole growth. 
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Eddington limit remains the relevant limit until the optical 
depth to electron scattering becomes so large that the pho- 
ton diffusion time is larger than the gas inflow time, in which 
case the photons are advected into the black hole along with 
the gas in a radiatively inefficient accretion flow. 

However, the situation is quite different when the total 
opacity is dominated by absorption. For UV photons emitted 
by the black hole, the force exerted on a shell of gas is: 



Figure 6. Cumulative time where the ratio of the black hole 
accretion rate to the Eddington rate is above the given value for 
the A2 simulation with different assumptions about the dust to 
gas ratio. Radiative forces on dust grains apparently prevent the 
black hole from reaching Eddington ratios greater than about 
20%. 
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/ 



Lmin(r, 1) 



(39) 



The incident UV photons deliver their momentum to the 
shell, after which they are converted to thermal IR photons 
where the opacity is much smaller. The thermal IR photons 
then simply leave the galaxy without further scattering un- 
til the galaxy becomes optically thick in the IR. The crucial 
point is that UV photons cannot build up and then diffuse 
out of the galaxy, as is the case when considering the electron 
scattering opacity and the dust opacity in the IR. The max- 
imum force that can be exerted on a shell of gas is f-[jvL/c 
where /uv is the fraction of the black hole bolometric lumi- 
nosity that is emitted in the UV. 

If the infalling clouds of gas are optically thin in the 
UV, then the force on a gas cloud is 



/ = 



LifvYTVY + TEs) Lt^s fuVt^VV 



«:es 



(40) 



where L is the bolometric luminosity and /uv is the fraction 
of the bolometric luminosity emerging in the UV. The cloud 
will feel a force much larger than if the dust opacity were 
not considered. However, if the cloud is optically thick in 
the UV (but not in the IR) then the force is 



/ 



L{fuv + TEs) 



(41) 



In this case, the force is nearly independent of the mass of 
the cloud; adding mass to an already optically thick cloud 
increases the gravitational force but not the opposing radia- 
tive force. This fact makes it difficult for the dust opacity in 
the UV to have a dramatic effect on the gas falling toward 
the black hole. 

Figure [9] shows the distribution of optical depths for 
the AO and A2 simulations and variants. The typical optical 
depths are far below unity, although a small fraction of the 
time the simulations become optically thick, at least along 
some lines of sight. However, the distribution of minimum 
optical depths for each simulation snapshot show that there 
are almost always optically thin lines of sight available that 
will allow IR photons to escape without experiencing multi- 



ple scatterings that are necessary for eq. ( 38 1 to be relevant 



In spite of dust having at best a minor effect on the dy- 
namics of the ISM, it has a dramatic effect on the radiation 
emerging from the galaxy. With no dust grain destruction, 
the optical depth in the UV from the center of the galaxy 
is typically of order unity. This means that the typical UV 
photon from the black hole is absorbed on its way out of the 
galaxy. Therefore observations of the AGN will be dramat- 
ically different for different assumptions about the creation 
and destruction of dust, even if the dust does not dramati- 
cally affect the dynamics of the galaxy. 

Absorbed UV/optical photons from the AGN are con- 
verted to IR photons, but the galaxy is typically quite opti- 
cally thin in the IR, so these photons escape without much 
fanfare. If the ISM were optically thick to IR photons, then 
they could have a dramatic effect as they scatter out of the 



galaxy, as pointed out by Murray et al. ( 2005 1 and Thomp- 



son et al. ( 2005 1, but we do not find the conditions necessary 



for this effect to prevail in our simulations. We find that in 
our detailed treatment of the radiative transfer, the factor 



r in eq. (38 1 never exceeds unity by a significant amount. 



cause the star light is emitted at larger radius and need 
not traverse the densest part of the ISM. The typical stel- 
lar UV photon is produced at roughly the half-light radius 
of the galaxy, from which point the typical optical depth is 
far less than unity. Therefore the stellar emission is largely 
unaffected by dust in the galaxy. Dramatically different as- 
sumptions about the dust to gas ratio have a minor effect 
on both the dynamics due to and the observations of stellar 
photons. 

A subtlety arises when considering the distributed ra- 
diation source (the stars) in the case where the ISM is opti- 
cally thick to absorption and the gas density is falling steeply 
(faster than r~^). In this case it is possible for the luminos- 
ity at a given radius to be negative, i.e. the net flux points 
inwards. If there were no absorption, then every ingoing pho- 
ton would eventually become an outgoing photon, either by 
traveling unimpeded and emerging at a different point on 
the same sphere of radius r, or else by scattering until it 
diffused to radius r. This must be the case or else photons 
would build up at radii smaller than r. 

However, when absorption is important, UV photons 
may propagate inward, undergo absorption, and re-emerge 
as IR photons. If the density profile is steep enough, this 
situation robustly leads to negative UV luminosities for radii 
smaller than the radius from which the typical UV photon is 
emitted. That is to say that the stars at larger radius exert 
a force acting to compress gas near the center of the galaxy. 
The force is not balanced by outgoing IR photons (which 
must after all carry the energy of UV photons absorbed at 
small radius) because the opacity is smaller in the IR. The 
momentum carried by the IR photons is deposited at larger 
radius. 

We find that this is often the case in our simulations: 
the gas density profile is generally steep enough to permit 
negative luminosities to develop within a few kiloparsecs of 
the center of the galaxy. This effect can be seen in Figures 
[7]and[8l 

There are situations in which dust plays a major role, 
such as galaxies undergoing starbursts or merger-induced 



quasar activity with copious quantities of cold gas ( Murray 



Most stellar UV/optical photons are not absorbed be- 



|et al.|2005| [Thompson et al.|2005|[Debuhr et al.|20l6||2011 
We have performed one simulation designed to be similar to 
an ultra-luminous infra-red galaxy having a star formation 
rate of several hundred solar masses per year. The simulation 
uses the A2 feedback parameters but starts with 3 x W^"Mq 
of gas, corresponding to a gas fraction of 10%, as could be 
induced by a major merger or an inflow of cold gas. The star 
formation rate reaches peaks of 200 solar masses per year, 
approaching the ULIRG conditions. 

Figure [lO] shows the distribution of Eddington ratios 
for different assumptions about the dust to gas ratio in this 
simulation. The simulation that includes dust creation and 
destruction is quite similar to the maximum dust case be- 
cause the short dust creation time ensures that most of the 
gas has a dust to gas ratio near that of the Milky Way, as 
is observed in the ULIRG case. Assuming that there is no 
dust in the simulation yields a dramatically different distri- 
bution of Eddington ratios, reaching much larger values. The 
time averaged black hole accretion rates are 0.045 Mq yr~^ 
for the case with maximal dust, 0.061 Mq yr~^ for the case 
with dust creation and destruction, and 0.32 Mq yr~^ for the 
case with no dust. Thus the presence of dust may suppress 
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Figure 9. In the upper two panels, the cumulative distribution of optical depths due to dust in the IR using the tracer field that tracks 
dust creation and destruction in the simulation to calculate the dust opacity. The lower two panels show the same quantity assuming 
a constant dust-to-gas ratio (corresponding to values appropriate for the Milky Way) so that the optical depth scales directly with the 
column density of the gas. For solid lines, each line of sight in each snapshot contributes separately to the distribution. Dashed lines show 
the distributions where only the line of sight with the highest or lowest optical depth from each snapshot contribute. The simulations 
only rarely become optically thick in the IR, and when this happens, the optical depth is only a few times unity. Furthermore, the 
distribution of minimum optical depths shows that there are always optically thin lines of sight available that allow photons to exit the 
galaxy without scattering multiple times. For the simulations labeled "No Dust," the simulation neglects any radiative forces associated 
with dust, but it nevertheless tracks the dust tracer field, allowing those simulations to be included in these plots. 



black hole growth by a factor of five to seven depending on 
the assumptions surrounding dust creation and destruction. 

The reason for this is that the conditions now resem- 



ble those envisioned by Thompson et al. (2005|) and Murray 



et al. ( 2005 1: the gas surface density is high enough that the 



galaxy is optically thick in the infrared. The infrared pho- 
tons build up in the galaxy until the gradient of the photon 
number density is sufBcient to carry the black hole luminos- 
ity by photon diffusion. In this case, the relevant opacity to 
limit accretion is not the electron scattering opacity but the 
dust opacity in the infrared. Dust then has a dramatic effect 
on the gas dynamics and on the growth of the black hole. 

Figure [TT] shows the cumulative distribution of optical 
depths due to dust in the IR for the ULIRG-like simulation. 
Along some lines of sight, the situation resembles that envi- 
sioned by Thompson et al. ( 2005 1 and Murray et al. ( 2005 1 



a significant but sub-dominant fraction of the time — optical 
depths are a few times unity and can exceed ten. However, 
the distribution of minimum optical depths for each simula- 
tion snapshot shows that IR photons would always be able 
to escape via optically thin "windows" . In order for equation 
|38| to pertain, the IR photons must be effectively trapped in 
the galaxy so that the only way for them to exit the galaxy 



is via diffusion. If there are optically thin lines of sight, this 
will not occur. 



6 DISCUSSION AND CONCLUSIONS 

We have outlined a method for computing the forces on in- 
terstellar gas due to absorption of photons by dust grains. 
The radiative transfer equation is recast into a set of dif- 
ferential equations by taking moments in the photon propa- 
gation direction. The set of differential equations comprises 
a boundary value problem, the solution of which gives the 
photon field which can be used to compute forces on the 
ISM. We identified a numerical method that reliably gives 
the solution to the set of differential equations even for very 
large optical depths. 

The method gives the correct asymptotic behavior in all 
of the relevant limits (dominated by the central point source; 
dominated by the distributed isotropic source; optically thin; 
optically thick to UV/optical; optically thick to IR; photon 
field nearly isotropic; photon field highly directed) and rea- 
sonably interpolates between the limits when necessary. The 
method is explicitly energy conserving so that UV/optical 
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Figure 10. Cumulative time where the ratio of the black hole 
accretion rate to the Eddington rate is above the given value 
for a simulation designed to resemble a ULIRG with different 
assumptions about the dust to gas ratio. For large enough gas 
column densities, radiative effects due to dust opacity have a large 
effect, and a careful treatment of dust production and destruction 
produces results similar to the "Max Dust" case, with "ULIRG" 
SMBH luminosity in the quasar range for several percent of the 
time. 




Figure 11. Cumulative distribution of optical depths due to dust 
in the IR for the simulation designed to resemble a ULIRG. The 
"Max Dust" case is shown, but high gas densities resulted in 
sufRciently fast dust creation that the other cases do not differ 
significantly. Typical optical depths along a typical line of sight 
are of order unity. There are almost always optically thick lines 
of sight where the optical depth is larger than ten a significant 
fraction of the time. However, even in this case there are always 
optically thin lines of sight that would allow photons to easily 
escape from the galaxy. 



photons that are absorbed are not lost, but are rather redis- 
tributed to the IR where they may scatter out of the galaxy. 
It does not require spherical symmetry, but it only handles 
the radial component of radiative forces. 

We have also described a physically motivated simpli- 
fication of the above method based on splitting the pho- 
ton field into ingoing and outgoing streams. This method is 
much simpler to implement and shares many of the desirable 
characteristics of the more complicated method. 

We have implemented both of these numerical meth- 
ods in a hydrodynamic code with source terms for energy, 
momentum, and mass appropriate to study black hole fuel- 
ing and star formation in early type galaxies. Both radiative 
and mechanical feedback are included for both the AGN and 
the stars. We allow for dust destruction by sputtering and 
dust creation by supernovae, planetary winds, and option- 
ally gas-phase grain growth. 

In the case of secular fueling of a normal SMBH in an 
elliptical galaxy, the dynamics of the simulations are not 
greatly affected by the presence or absence of dust. At high 
AGN mechanical feedback efficiencies, dust has very little 
effect because mechanical feedback dominates. At lower me- 
chanical feedback efficiencies, dust has some effect but large 
changes in our assumption about the dust-to-gas ratio lead 
to measurable but relatively modest changes in the black 
hole growth, star formation, and galactic winds. This is true 
even for extreme assumptions about dust, ranging from no 
dust at all to assuming dust-to-gas rations appropriate to 
the Milky Way pertain at all times. Allowing for the possi- 
bility of dust destruction by sputtering in hot gas makes the 
simulations behave similarly to the "No Dust" case. 

If the galaxy is optically thick to IR photons, repeated 
scattering of these photons can have a dramatic effect on 
the momentum balance of the ISM. We have verified this 
by simulating a galaxy with a very large mass in the ISM, 
as might follow a major merger with a gas rich galaxy, in 
which case the presence or absence of dust makes a large 
difference, as expected. 

However, we do not find that this situation is typically 
realized in our simulations. The opacity of dust in the IR 
is about five times the electron scattering opacity, so a sys- 
tem must be Compton thick before this scattering starts to 
have a large effect. Our simulated galaxies are occasionally 
Compton thick, but no more, and they typically have much 
lower column densities. Even with no dust destruction, we 
find that the optical depth in the IR is at most a few times 
unity, and it reaches that value only for a small fraction of 
the total simulation time. We find that it is very rare for 
the optical depth in the IR to significantly exceed unity, as 



is assumed by Debuhr et al. (2010 2011 1. When dust de- 



struction is allowed, the optical depths in the IR is typically 
smaller by several orders of magnitude. 
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APPENDIX A: THIRD ORDER RADIATION 
TRANSPORT 

The quality of a solution to the radiation transport equa- 
tion obtained by taking moments depends on the accuracy 
of the physical assumptions implicit in the chosen analytic 
form for the specific intensity. If those assumptions are not 
satisfied, taking the expansion to higher order is unlikely 
to provide benefits commensurate with the costs in terms 
of complexity and computational resources. We believe that 
the assumptions embodied by eqs. pi and ( 12 1 are quite rea- 
sonable for the problem at hand. Nevertheless, it is possible 
to avoid using eq. ( |12[ ) by taking one additional moment of 
the specific intensity and radiation transport equation. We 
provide the equations here for reference but do not pursue 
the method further in the present paper. 

Define the third moment of the specific intensity to be: 



Q 



^ /(/i, (;/))dp. 



Given eq. ([3| we find 

^ 47rB 
Q = -—- + TvD . 
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intensity and four moments of the mean intensity, so we can 
solve for the closure relation: 

_ 47r J + SF + 3cP 
Q- g . 

To obtain the differential equations, we take the equa- 
tions for the zeroth and first moment derived above, supple- 
ment them with the second moment of the radiation trans- 
port equation, and then insert the closure relation to obtain 
the following set of three equations: 
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Adopting the use of these equations and closure rela- 
tion is perhaps formally more appealing than the method 
outlined in Section [S] However, the cost is one additional 
equation to integrate and perhaps less obvious physical con- 
tent to the closure relation and the equations themselves. 
We have included the present equations for reference, but 
do not use this method of obtaining solutions to the radia- 
tion transport problem in the present work. 



APPENDIX B: 
TRANSPORT 



SIMPLIFIED RADIATION 



In Sectionlslwe have described a complete algorithm to solve 
for the photon field and, subsequently, the radiative forces 
on interstellar gas. The algorithm is exact in the sense that 
if the physical approximations used in its derivation are cor- 
rect, the method will compute the correct photon field and 
radiative forces. 

However, solving for the radiation field dominates the 
computational resources required for the calculation pri- 
marily because the system differential describes a bound- 
ary value problem rather than an initial value problem and 
hence finding the solution requires iteration. We would like 
a another solution method that is computationally efficient 
while respecting the essential physics of the previously de- 
fined algorithm. 

Consider the problem where the only source of radi- 
ation is the black hole and the scattering opacity is zero. 
The UV/optical photons are always radially outgoing at all 
times. Some UV/optical photons are converted to IR pho- 
tons which diffuse out of the gala^cy, but the UV/optical 
radiation field never becomes isotropic. In this case the fiux 
is simply related to the mean intensity {F = An J) and we 
can eliminate eq. (Ill to obtain the single equation for the 



UV and optical components 

dL 

dr 



-pHaL . 



(Bl) 



We have three constants in the assumed form of the mean 



The boundary condition is L{r min) ~ Lbh- This is now an 
initial value problem rather than a boundary value problem, 
and can be trivially solved by integrating from the inner edge 
of the simulation grid. 
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This change greatly simphfies the solution procedure. It 
would be very useful to obtain a simple approximate treat- 
ment of this problem that as an initial value problem rather 
than a boundary value problem. We seek such an algorithm 
presently. 



by solving the boundary value problem for one line of sight 
in one simulation snapshot. The total optical depth along 
the line of sight is 13.5. The solution given by the simplified 
scheme is not perfect, but it is remarkably good given the 
large optical depth and the reduced algorithmic and com- 
putational cost. 



Bl Differential equations 

Conceptually, there are two streams of radiation: an ingoing 
stream and an outgoing stream. If the optical depth is low, 
then the photons in the ingoing stream are likely to success- 
fully traverse the inner parts of the galaxy and emerge as 
outgoing photons. In this case, all of the radiation emitted 
by stars should be considered to be added to the outgoing 
stream. Doing so gives the correct forces on the material in- 
side and outside the emitting sphere due to the well-known 
result that for a inverse square force law, there is no net 
force on a test particle inside a spherical shell of material; 
and a test particle outside the spherical shell feels a force 
equal to that of a point source of equivalent mass, charge, 
or luminosity, no matter whether the force is gravitational, 
electrostatic, or due to momentum transfer by photons. 

If the optical depth is large, then the ingoing photons 
are likely to be absorbed in the inner part of the galaxy. 
Then only half of the emitted photons should be added to 
the outgoing stream. The other half should be added to the 
ingoing stream, where they will in due course be absorbed. 

This leads to the following two equations: 

dLout ^ 2 ,T, 7^ r 

= 4nr wE - pnaLout , 



dr 
and 
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dr 
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where ^ is the fraction of photons emitted at a given radius 
that are likely to become outgoing photons, either because 
they were emitted in the outward direction, or because they 
were emitted inward but are likely to successfully traverse 
that region without being absorbed. 



B2 Solving the differential equations 

The boundary conditions are I/out(fmin) ~ and 
iin(''max) = 0, SO that eq. (|B2[) is solved by integrating out- 



ward and eq. ( B3 1 is solved by integrating inward. 



Summing eqs. (B2| and (B3l simply gives eq. (10 1, as 



expected. The advantage of splitting the radiation into an in- 
going stream and an outgoing stream is that it transforms a 
boundary value problem into an initial value problem, lead- 
ing to a great reduction of the effort required to find a nu- 
merical solution. The quality of the solution will depend on 
the accuracy of the function ^P is estimating the contribu- 
tion that the stars at a given radius make to the ingoing 
and outgoing streams of radiation. This can of course be 
checked by comparing to the solution obtained by solving 
the full boundary value problem. 

We have found that this scheme gives reasonable results 
at a fraction of the algorithmic complexity and computa- 
tional cost of obtaining the exact solution. Figure [BT] shows 
the results of solving the radiation transfer equation by the 
simplified two-stream method and the full solution obtained 



B3 Concerning ^ 

If Pabs is the probability that a given ingoing photon is ab- 
sorbed before again reaching the radius at which it was 
emitted (and hence becoming an outgoing photon) and 
Ptrans = 1 — Pabs is the probability that a given ingoing 
photon is transmitted so that it again arrives at the radius 
at which it was emitted and becomes an outgoing photon, 
then: 

_ Pabs _ 1 + Ptrans 
2 ^ 2 



* 



(B4) 



The probabilities for absorption and transmission appear 



in eqs. ( B2 1 and ( B3 1 only through the quantities '^ and 
1 — ^. Note that of the four ways of writing ^ and 1 — VE" in 
terms of transmission and absorption probabilities, pabs and 
Ptrans are added to an order-unity constant in all cases ex- 
cept one: 1 — 'i' — pabs/2. An argument could be made that 
one should take special care to obtain the correct asymptotic 
behavior for Pabs when pabs <C 1 in order to avoid large frac- 
tional errors in the ingoing radiation stream. However, this 
is the optically thin limit, exactly where the ingoing radia- 
tion stream is expected to be of limited importance. There- 
fore any expression for ^ with the correct asymptotic limits 
and a transition between the two limits when the system 
goes from optically thin to optically thick conditions should 



introduce only limited fractional errors into eqs. ( B2 \ and 



(B3l 



The simplest estimate of the fraction of photons emitted 
at a given radius that are likely to be absorbed at smaller 
radii comes from finding the apparent size of the sphere 
where r = 1 as seen from the radius r, where the photons 
are emitted. 



* = 1 



1 + e- 



rl 



max(rf. 



(B5) 



where r is the optical depth from r to infinity, t = J pn dr, 
and ri is the radius where the optical depth is unity: 1 — 
J pndr. Here the factor 1/(1 + e""^) interpolates between 
the correct asymptotic expressions for the case where r ^ 
1 , ri ^ r, in which $ ~ 1 — O.Sr^/r^ and the case where 
r > 1 , n > r in which * ~ 1/2. 

More rigorous formulae can be derived assuming spher- 
ical symmetry and a power- law dependence of the mean free 
path for absorption on the radius. Given the geometry shown 
in Figure |B2[ the optical depth as a function of the angle a 
between the direction of the ray and the line connecting the 
emission point to the center of the coordinate system is 



T{a) = 



r sin a dp 






:+/J) 



(B6) 



where r is the radius of emission, l{s) = 1/ p{r)K,{r) is the 
mean free path to absorption and Q, /?, and s are defined by 
Figure [B2] The formula gives the optical depth for a photon 
to traverse the chord of the circle in Figure |B2] 
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Figure Bl. The luminosity as a function of radius for the full relaxation based algorithm (Section |3.2| and the simpler two-stream 
based algorithm (AppendixlBj. Dashed lines indicate inward-directed luminosities. The solution given by the two-stream approximation 
described in Appendix^ has about 50% more flux at infinity than the "true" solution, and it misses the true solution around 100 pc. 
However, the solution is quite good given the gains in algorithmic simplicity and reduced computational costs. The function ^ associated 
with C = 2 was used for this figure, where Q is the index of the assumed power law dependence of photon mean free path on radius used 
in defining ^ (see Appendix |B3| for details). 




Figure B2. Geometry for calculation of optical depth in Section 



If the mean free path as a function of radius can be 
approximated by a power law 



i{s) ^lA- 



«(?)' 



C>o 



(B7) 



where Iq is the local mean free path at the radius r at which 



the photons are emitted, then eq. ( B6 1 can be written in 



terms of hypergeometric functions. However, assuming that 
C, is an integer gives the simple formulae involving elemen- 
tary functions given in Table [BT] 

The transmission probability ptrans is then 



r/2 
Ptrans = / siu Of fi" 

'0 



(B8) 



Table Bl. Optical depth as a function of photon emission direc- 
tion in the case of a power-law mean free path, where C, is the 
power-law index. 



C <a) 



(2r//o) coso 

(2r/io)logcot(Q:/2) 

{r/l(}){TT — 2a)/sina 

(2r//o) cos a/ sin^ a 

(r/2lo){-K — 2a + sin2a)/sin^ a 



so that ^f can be recovered from eq. ( B4 1 



The integral in eq. ( B8 1 unfortunately cannot be ex- 



pressed in terms of elementary functions. However, for a 
given dependence of mean free path on radius, it is a func- 
tion only of the radius in units of the local mean free path 
r/k). For integral values of (^ the integrand consists of ele- 
mentary functions, so it can be tabulated at the beginning of 
a calculation and used throughout with little computational 
cost. For values of ^ from to 4 (corresponding to cases 
from p = constant to p oc r~* if the opacity is constant), the 
values of ptrans and pabs do not depend sensitively on (,. 

Figure |B3| shows pabs and ptrans as the as a function of 
optical depth for different values of ^. The precise value of 
^ makes a moderate difference when the mean free path is 
roughly equal to the radius of the sphere. In the optically 
thick limit, the asymptotic form of ptrans is independent of 
(" and is given by 



PtT 



0.005 



/ ^0 

UOO 



for /() < 7" . 



(B9) 
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For the optically thin case, the asymptotic value of Pabs de- 
pends on ^. For values of (^ from to 4, the asymptotic form 
can be approximated as: 



Pabs 



1 

Too 



1 + 6It 



(D ( 



'l+0.4(C/4)'^ 



for /o>r-.(B10) 



Following the dictum discussed above that errors in 
Ptrans are tolerable because ptrans is always added to an or- 
der unity constant, while the same is not true for Pabs for 
the ingoing radiation stream when pabs is small, then one 
may take the approximation: 



Pabs 



Pabs 
1 + Pabs 



(Bll) 



with Pabs given by eq. ( BIO I over the whole range of optical 



depths. The value of 'I/ may then be recovered using equation 
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Figure B3. Probability of transmission and absorption for spherically symmetric systems with power-law mean free paths of index ^. 
Solid lines give Pabsi dashed lines give ptrans- 
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